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The thickness of Hi in galactic discs under MOND: theory 
and application to the Galaxy 
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ABSTRACT 

The outskirts of galaxies are a very good laboratory for testing the nature of the 
gravitational field at low accelerations. By assuming that the neutral hydrogen gas 
is in hydrostatic equilibrium in the gravitational potential of the host galaxy, the 
observed flaring of the gas layer can be used to test modified gravities. For the first 
time we construct a simple framework to derive the scaleheight of the neutral hydrogen 
gas disc in the MOND scenario and apply this to the Milky Way. It is shown that using 
a constant gas velocity dispersion of ~ 9 km s -1 , MOND is able to give a very good 
fit to the observed H I flaring beyond a galactocentric distance of 17 kpc up to the last 
measured point (~ 40 kpc). Between 10 and 16 kpc, however, the observed scaleheight 
is about 40% more than what MOND predicts for the standard interpolating function 
and 70% for the form suggested by Famaey & Binney. Given the uncertainties in the 
non-thermal pressure support by cosmic rays and magnetic fields, MOND seems to be 
a plausible alternative to dark matter in explaining the Milky Way flaring. Studying 
the flaring of extended H I discs in external edge-on galaxies may be a promising 
approach to assess the viability of MOND. 

Key words: galaxies: haloes galaxies: kinematics and dynamics - galaxies: structure 
- gravitation - Galaxy: disc 



1 INTRODUCTION 

The success of Modified Newtonian Dynamics (MOND) in 
reproducing the observed velocity rotation curves of spiral 
galaxies without dark matter haloes is amazing. Only in 
about 10% of the roughly 100 galaxies considered in the 
context of MOND does the predicted rotation curve differ 
significantly from that observed (e.g., Sanders & McGaugh 
2002). Bekenstein (2004) has suggested a theory that leave a 
door open to embed MOND within a physically acceptable 
dynamical framework. However, other aspects of MOND in 
spiral and elliptical galaxies are yet to be investigated. It 
would be desirable to see if MOND is able to reproduce the 
full three-dimensional structure of galaxies (e.g., Kuijken & 
Gilmore 1987; Hernquist & Quinn 1987; Buote & Canizares 
1994; McGaugh & de Blok 1998; Milgrom 2001; Stubbs & 
Garg 2005; Nipoti et al. 2007). 

The flaring of the H I gas layers has been used to deter- 
mine the Battenings of dark haloes for a few external galax- 
ies because it is sensitive to the shape of the dark matter 
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halo (e.g., Oiling 1996; Becquaert & Combes 1997). The al- 
ternative theories to dark matter may also have observable 
consequences on the vertical structure of galactic discs. For 
instance, studies of the flaring of the H I discs have been able 
to rule out the magnetic alternative to dark matter because 
the magnetic pressure would produce extraordinary thick 
discs (Cuddeford & Binney 1993; Sanchez- Salcedo & Reyes- 
Ruiz 2004). In modified gravities such as MOND, where the 
gravity comes mainly from the disc, the potential may be 
much flatter leading to far more thin discs (e.g., Sanchez- 
Salcedo & Lora 2005). In this paper we suggest the usage of 
the flaring of Hi discs to test modified gravities. 

The Milky Way is an attractive target to carry out self- 
consistent tests because of the excellent flaring data, large 
radial extend and rather well-known distribution of stars and 
gas. In the conventional Newtonian view, Kalberla (2003) 
and Kalberla et al. (2007) have found no easy way to match 
a standard flattened dark halo to the Galactic flaring data. 
They discuss the role of support by non-thermal pressures, 
non-axisymmetric mass distributions of dark and baryonic 
matter, spiral structures and so on, and conclude that, be- 
sides a massive round dark halo, a self-gravitating dark mat- 
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ter disc with mass ~2x 10 Mq plus a dark matter ring at 
13 < R < 18.5 kpc with mass ~ 2.5 x 10 10 M Q are required. 
This result is especially intriguing given other estimates for 
the shape of the Milky Way's halo that suggest that it could 
indeed be rather spherical or even prolate (e.g., Helmi 2004; 
Johnston et al. 2005; Law et al. 2005; Belokurov et al. 2006; 
Fellhauer et al. 2006). In addition, the existence of a mas- 
sive dark matter disc is unexpected in the collisionless CDM 
paradigm. This might indicate that more physics should be 
added in the simulations of structure formation (e.g., bary- 
onic physics) or the lack of complete modelling of the flaring 
problem. 

It is commonly argued that MOND is a more falsifi- 
able theory than CDM hypothesis because it has less ad- 
justable parameters to fit astronomical observations. One 
expects that if it is not easy to explain the flaring in CDM 
theory, it should be even more difficult in MOND. It is in- 
teresting to see if, given the measured baryon distribution 
of the Milky Way, MOND is able to successfully explain the 
Galactic flaring. 

In section 2, we discuss the vertical distribution of self- 
gravitating mass layers under MOND. In particular, we give 
the semi-analytic framework to derive the radial dependence 
of the thickness of the Hi layer in axisymmetric galactic 
discs. In section 3, we apply our method to derive the ex- 
pected MOND Hi flaring in the Milky Way. By comparing 
the predictions with recent determinations of the thickness 
of the H I disc in the Galaxy, we are able to test the ability 
of MOND to reproduce the H I flaring in an extended radial 
range (10 kpc < R < 40 kpc). 



2 FLARING OF DISCS: THEORETICAL 
APPROACH 

2.1 Background and motivation 

A possible interpretation of the missing mass in spiral galax- 
ies is that the flatness of the rotation curves reflects a failure 
of the Newton law rather than the existence of vast quanti- 
ties of dark matter. Milgrom (1983) suggested that for ax- 
isymmetrical systems, the real acceleration at the midplane 
of the configuration, g, is related with the Newtonian accel- 
eration, Qn, by: 



MI<?l/ a o)<7 = Sat, 



(1) 



where ao is an universal acceleration of the order of ~ 10~ 8 
cm s~ 2 , and p(x) is a monotonic and continuous function 
with the property that n(x) ~ x for i < 1 (deep MOND 
regime) and p(x) ~ 1 for x S> 1 (Newtonian regime). Al- 
though some effort has been made to derive the form of the 
interpolating function from basic principles (e.g., Bekenstein 
2004), it is arbitrary but once it is fixed from astronomical 
observations, it is universal. A rather general form is 



p(x) 



(i + x n yi n ' 



(2) 



where n(x) > and obeys that n'(x) < for x <§C 1 and 
n'(x) > for x ^> 1. The simplest case is to assume n a pos- 
itive real constant. The exponent n determines how abrupt is 
the transition between the Newtonian to the MOND regime. 
Since for n > 2 the interpolating functions do not differ 



much, Milgrom (1983) selected n = 2 for simplicity. Interest- 
ingly, the standard (n = 2) interpolating function succeded 
in fitting the rotation curves of most external galaxies, be- 
ing the stellar mass-to-light ratio of the disc the only free 
parameter (e.g., Sanders 1996; Sanders & McGaugh 2002; 
Sanders 2004). Recently, Famaey & Binney (2005, hereafter 
FB) suggested that n = 1 is better for MOND to fit the 
inner rotation curve of the Basel model of the Milky Way 
(Bissantz et al. 2003). For a sample of galaxies having a 
gradual transition from the Newtonian limit in the inner 
regions to the MOND limit in the outer parts, Famaey et 
al. (2007) and Sanders & Noordermeer (2007) conclude that 
with respect to rotation curve fitting alone, there is no rea- 
son to prefer n — 1 to n = 2, it is only the plausibility of 
the relative M/L values for bulge and disc as well as the 
generally smaller global M/L, that support n = 1. More- 
over, the best fit for the sample of Famaey et al. (2007) was 
obtained for n = 1. From the \ 2 statistics of the global fits, 
they derived n > 0.85 at one-sigma uncertainty. Since values 
n > 2 seem to be incompatible with the measured terminal 
velocities interior to the Solar circle and require disc M/L 
ratios slightly higher than the predictions of stellar popu- 
lation synthesis models (McGaugh 2005), it turns out that 
values in the range 0.85 < n < 2 are allowed by current 
observations, with n = 1 the most favoured value. 

The phenomenological formulation of MOND by Eq. (1) 
cannot be applied to an arbitrary self-gravitating system. 
Bekenstein & Milgrom (1984) suggested a Lagrangian the- 
ory and also a nonlinear differential equation for the nonrel- 
ativistic gravitational potential produced by a mass density 
distribution p: 



V$, . .... , 

p I - — '- 1 v$ 

fflo 



= 4nGp. 



(3) 



Only for very special configurations (one-dimensional sym- 
metry -spherical, cylindrical or plane symmetric systems- 
or Kuzmin discs) the MOND field is related to the Newto- 
nian field by the algebraic relation Eq. (1) (Brada & Mil- 
grom 1995). The modified Poisson equation given by Eq. (3) 
overcomes many of the conceptual maladies of the algebraic 
equation but it is very difficult to solve. 

In a MONDian baryonic Universe p is the density dis- 
tribution of the seen baryonic matter. A Newtonist, who 
assumes that $ solves the Poisson equation, will deduce a 
fictitious dark matter density: 

p dm = (4ttG)- 1 V 2 $-p. (4) 

As Milgrom (2001) demostrated, MOND predicts for disc 
galaxies a distribution of fictitious dark matter that com- 
prises a dark disc and a rounder halo. The surface density 
of fictitious dark matter in the disc is 



p(g+/a ) 



I E, 



(5) 



where E is the real (baryonic) surface density of the disc 
and g + is the total MOND acceleration just outside the disc. 
At large radii where p <C 1, the 'dark disc' dominates, im- 
plying that Newtonian disc-flaring analyses should indicate 
highly flattened oblate dark haloes. Although the initial ap- 
plications of the flaring technique for NGC 891 and NGC 
4244 implied highly flattened dark haloes, the gas layer flar- 
ing method does not return systematically highly flattened 
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haloes (Oiling & Merrifield 2000; Merrifield 2002; Narayan 
et al. 2005 for the case of the Milky Way). The question 
that arises is: are the predictions of MOND compatible with 
measurements of the flattening of galaxies? Since it is not 
possible for MOND to mimic prolate dark matter haloes, 
a clear-cut probe of the existence of a prolate halo in an 
isolated galaxy would be difficult to reconcile with MOND. 



2.2 Scaleheight of one-dimensional planar mass 
layer in deep MOND 

Assuming pure planar symmetry, the only component dif- 
ferent from zero is g z . The vertical hydrostatic equilibrium 
equation is: 



dP d$ 
-P- 



dz 



dz 



(6) 



Integrating this equation from to z and exploiting MOND 
field equation (Eq. 3) , we find 



In the deep MOND, p ~ a 1 d<fr/dz, and hence the RHS of 
the above equation can be integrated immediately: 



P(z) - P(0) 



dl> 



6irGao V dz 



(8) 



The value of the vertical force d<fr/dz can be inferred again 
from the MOND field equation after integrating it between 
±z, d$/dz = \/27rGaoE(z), where E(z) is the integrated 
surface density £(z) = f* z pdz'. Substituting this expres- 
sion into Eq. (8), 



P(z) - P(0) = -- (2^Ga„) 1/2 E(z) 3/2 . 



(9) 



For an isotropic and isothermal gas with one-dimensional 
velocity dispersion a, the midplane pressure is P(0) = poo , 
where po is the density at z = 0. Adopting the boundary 
condition P{oo) — 0, we derive the pressure at the midplane: 



poa 2 = i(2^Ga ) 1/2 E 3 / 2 , 



(10) 



with Eoo = E(oo). When the scaleheight zo is defined as 
zo = Eoo/(2po), it satisfies 



zo = 



2 V27rGa E o 



(11) 



for our isothermal self-gravitating layer in the deep MOND 
regime. The dependence of zo on Eoo is obviously different 
than in the Newtonian case. In order to facilitate comparison 
it is convenient to rewrite this formula in a more familiar way 
exploiting the relation p x = p(oo) = <Xq 1 v / 27rGaoE 00 , 



zo = 



3 Hoo<y 
IttGE^ 



(12) 



Therefore, the usual recipe of replacing G — > G/p to obtain 
the MOND value from the Newtonian expression overesti- 
mates the scaleheight by a factor 4/3 in this case. 



2.3 Infinite mass sheet in an external 
gravitational field 

Suppose now that our mass layer is under a constant exter- 
nal field g ext = (g x , 0, 0) and that the dynamics is dominated 
by g x , i.e. g x 3> g z , where g z is the vertical acceleration cre- 
ated by the self-gravity of the layer. In that case the modified 
Poisson equation (Eq. 3) reads 



p= (4ttG)" 



dz V ^ dz 



(13) 



with p,(g/ao) and g = (gi + gl) 1 ^ 2 ~ gx in an external- 
dominated field. Therefore p is constant (independent of x 
and z), and thus 



P = 



P 



d 2 $ 



4ttG dz 2 



(14) 



Following the same steps as in the previous subsection, the 
vertical scaleheight of the layer in an external dominated 
field is 



zo = 



pa 



7rGE 



(15) 



Note that p(g x /ao) =constant. When g x 3> do, we have 
p w 1 and we recover Spitzer's formula. If, on the opposite 
case, the layer is in the deep MOND regime (g z <C g x <S cio), 
then p ~ g x /ao and zo is proportional to g x . This result has 
interesting consequences in galactic discs. Consider now a re- 
gion of the disc outside the central region such that gn^> g z . 
If ga is mostly determined by the inner stellar disc and g z 
by the local self-gravity of the H I layer, then gn plays the 
role of the external field. According to the ongoing discus- 
sion, if the adopted stellar mass-to-light ratio is increased 
and thus its mass, then the H I layer in that region becomes 
thicker. This result is reminiscent of the external field effect 
(Milgrom 1995) and will be relevant in §2.5. 



2.4 Gravitational potential in axisymmetric thin 
discs 

The MONDian Poisson equation (3) is difficult to solve even 
numerically. However, it can be simplified when studying the 
vertical structure of the outer parts of axisymmetric galactic 
discs. The field equation can be written as: 



p(a:)V 2 <I> + Vp(x) ■ V<3> = 4ttGp, 



(16) 



where x = |V<E>|/a . For an axisymmetric configuration 
Equation (16) can be expanded in cylindrical coordinates 
(R, z) as: 

^)V 2 * + (||||)(l + e)=4 7r Gp, (17) 
where 

c = e*e#. (is) 

with 

«*•>-(£)(£)"'-£• 

and 

«*■>-(£) (£)"'- (£)(&)"'• ^ 
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Note that in a smooth disc galaxy gn < 0, and dg/dR < 
(see the Appendix), whereas sign((/ z ) = — sign(dgfdz) = 
— sign(z), implying £ < 0. No approximations at all have 
been made so far, and Eq. (17) is fully equivalent to the 
original one for axisymmetric systems. 



The factor ff^fjf in Eq. (17) can be rewritten as: 



d^_dp_ 
dR dR 



dp 9$ dx 
dxdR~d~R 



i_dpd^__d_ 

ao dx dR dR 



(21) 



In the last equality we have used that < 1 within H I 
galactic discs (see the Appendix), so that x = |V$|/ao = 
x{R,Q)y/l +£| ~ x(R,0) =Vc(R,0)/(Ra ). Using that 



1 d 



RdR \ OR 



R 



9$ 



+ 



a 2 $ 



1 dv 2 d 2 $ 
R ~d~R + Ite 2 ' 



(22) 



we are lead to the final modified Poisson equation under 
MOND dynamics: 



d 2 $ 



n i y2 c d ( y2 c\ a* dv l (2 o\ 

( U a RdxdR\R J RdR' { ' 



where v 2 = v 2 (R, 0) is the circular velocity at the midplane 
of the disc. Or, equivalently, 

= 471 



-^ = ^ p -^ + ^ L ^dR{ii)-R^m- 



(24) 



where L(x) = d In p/d In a; > is the logarithmic deriva- 
tive of p at the same value of the argument. For a nearly 
constant rotation curve, the third term of the RHS is very 
small. Combining Eq. (24) with the hydrostatic equilibrium 
equation (Eq. 6), one can obtain the vertical distribution for 
an isothermal disc. 



2.5 Interpretation of the modified Poisson 
equation 

In order to gain some insight, below we compare the different 
terms of Eq. (24) with the corresponding terms in the New- 
tonian equation. For an axisymmetric disc, the Newtonian 
Poisson equation is 



9 2 $ 1 dv 2 

d^ = ^ Gp -RdR- 



(25) 



The equation above has been used vastly to study the ver- 
tical equilibrium in the standard dark matter scenario (e.g., 
Oiling 1995, 1996; Oiling & Merrifield 2000; Narayan et 
al. 2005). Obviously, the Newtonian equation can be recov- 
ered from the MONDian equation taking the limit ao — > 0. 
In fact, in the Newtonian limit, p — » 1, dp/dx — > 0, hence 
L(x) -> and Eq. (24) simplifies to Eq. (25). 

Comparing Eqs (24) and (25), we see that in the MOND 
case, the surface density is replaced by the effective surface 
density S//u, as expected from our discussion in sections 2.1 
(see Eq. 5) and 2.3. The confining effect of this term is in- 
versely proportional to p. As anticipated in section 2.3, if we 
increase the stellar mass of the disc by adopting a bit larger 
stellar mass-to-light ratio, the rotation velocity v c will in- 
crease correspondingly and hence, the effective gravitational 
constant G/p oc G/« 2 becomes smaller at these radii. 

Interestingly, Eq. (24) contains a term proportional to 
L(x); we will refer to it as the L-term. In principle, this 
term has no analog in Newtonian dynamics. In order to un- 
derstand its nature, consider the region of a certain galaxy 



where the rotation curve is flat, i.e. v c (R) =const, and write 
Eq. (24) in terms of the vertical and angular frequencies: 



v 1 = 47T-P+ (1 + OL(x)Q 2 
A 4 



(26) 



where £1(R) = v c /R. The L-term becomes important in the 
vertical equilibrium configuration when 



{l + Z)VL(x)L(x)Sr>4*Gp. 



(27) 



Centainly, the L-term is crucial when the centripetal accel- 
eration is sufficiently low that we are in the MOND regime, 
H(x) w x, and for midplane densities lower than the critical 
density given by: 

^) = ( 1 + ^4§r ^ 

Since the asymptotic velocity in MOND is vt = GMao, 
where M is the total mass of the disc, the critical density 
can be rewritten as p c (R) = (1 + ^)M/(4tvR a ). 

As a reference number, the critical density for a galaxy 
with v c = 220 km s" 1 is 2.3(1+0 x 10~ 3 M Q pc~ 3 at R = 20 
kpc. For Galactic parameters, we find that |£| w 0.03 at 
20 kpc (using our estimate, Eq. (A6), with a HWHM scale 
h = 500 pc and a = 9 km s _1 ). Since the midplane Hi 
volume density at 20 kpc is ~ 0.9 x 10~ 3 M pc~ 3 , we 
expect the L-term to be important or even dominant at the 
very outer disc, say R > 20(-u c /220kms~ 1 ) kpc. At these 
radii, £ can be neglected as compared to 1 and Eq. (26) can 
be written as: 



z/ 2 = 47r- y^pi + n 2 , 



(29) 



where we have explicitly quoted that the baryonic density 
is the sum of three components, the stellar, atomic (He cor- 
rected) and molecular hydrogen mass densities. Stubbs & 
Garg (2005) have proposed to use jointly the rotation ve- 
locity and velocity dispersion of galactic discs to test the 
self-consistency of MOND. From our analysis, we learn that 
p can be derived from the observed vertical velocity disper- 
sion a z , the scaleheight of the disc h and the surface density 
S, according to the following relation 

a'i = T^Sft+ ^L(x)Q, 2 h 2 . (30) 

2.6 L-term and the dark matter halo 

Imagine a Newtonist and a MONDianist trying to explain 
the dynamics of a particular galaxy, i.e. the same Since 
they use different field equations, they infer different mass 
distributions. Newtonists add a spherical dark halo of the 
form 

„,2 



Ph(r 



(31) 



47rGr 2 

to explain the flatness of the rotation curves in the outer 
parts of disc galaxies. In the following we argue that the L- 
term is the MOND counterpart of the Newtonian spherical 
dark halo. We start noting that, with this dark halo, the 
conventional Poisson equation at large radii reads: 



a 2 $ 



— =47rG[J2 Pi + p h 



(32) 
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Substituting for p h into Eq. (32) with v 2 — Q. 2 r 2 we obtain 



dz 2 



or written in terms of the frequencies: 

3 

v 2 =47rG^pi + fi 2 . 



(33) 



(34) 



Comparing Eqs (29) and (34) we see that, in the outer parts 
of galaxies where the dynamics is in the deep MOND regime, 
the L-term mimics the vertical pull created by the spheri- 
cal dark halo in the traditional Newtonian scenario (see also 
Milgrom 2001). In summary, MOND predicts for disc galax- 
ies a fictitious round dark halo and a fictitious dark disc. 
The resultant flattened potential along the vertical axis im- 
plies that, if MOND is correct, Newtonists should only infer 
oblate dark matter haloes around disc galaxies. The infer- 
ence of prolate (q > 1) dark haloes in isolated spiral galaxies 
would suffix to rule out 'classical' MOND. Of course, MOND 
is not fundamentally incompatible with the existence of un- 
seen matter and one might argue that an unseen prolate 
spheroid of stars can account for the Newtonian inference of 
a prolate dark halo, but this would weaken its appeal. 



3 APPLICATION: THE MILKY WAY 

One advantage of studying our Galaxy as compared to other 
external galaxies is that we know quite well its baryonic 
mass. For no other galaxy do we have information of compa- 
rable quality. Using the terminal velocities, Famaey & Bin- 
ney (2005) were able to put constraints on the form of the 
interpolating function given the measured baryon distribu- 
tion. In addition, the new All Sky Hi Survey (the Leiden- 
Argentine-Bonn [LAB] survey as it is known; Kalberla et 
al. 2005) provides us a reliable three dimensional structure 
of the H I disc of the Milky Way up to a minimum of 40 kpc 
from the Galactic centre. Currently, this is the most sensitive 
Milky Way H I line survey with the most extensive coverage 
both spatially and kinematically. Since the stellar contribu- 
tion drops down rapidly in the outer region, this huge H I 
disc can serve as an excellent probe for testing alternative 
theories to dark matter and thereby the nature of gravity. 
In this analysis we use MOND (as opposed to dark matter) 
as a possible candidate for reproducing the thickness of the 
H I layer of our Galaxy and hence test its standpoint in this 
regard. 



high velocity dispersion (12 km s _1 ). Therefore, the assump- 
tion of constant velocity dispersion should be reasonable af- 
ter averaging over azimuth. Under these assumptions, the 
thickness of the Milky Way's gas layer is dictated by the 
hydrostatic balance between the pull of gravity toward the 
Galactic plane and the pressure forces acting on the gas: 



idfh 
dz 



-Pi 



<9£ 
dz' 



(35) 



If the interstellar gas is in the form of small dense clumps, 
the vertical component of its velocity dispersion should be 
smaller than the radial component, in the same manner than 
stars. Even so, the term due to the tilt of the velocity dis- 
persion tensor is expected to be very small (Oiling 1995; 
Becquaert & Combes 1997; Kalberla 2003). The gradient 
of the mean pressure of the hot gas (T > 10'' K) is neg- 
ligible within the thickness of the the Hi disc because the 
scaleheight of hot gas is much larger than the scaleheight 
of H I. For instance, the decrease of hot gas density away 
from the plane is consistent with a local scaleheight of 3.2 
kpc. Overpressured bubbles of hot gas formed by supernovae 
explosions may be a source of driving porosity in the ISM 
and turbulent motions in the Hi. Porosity may modify the 
hydrostatic properties of the gas layer. However, for a turbu- 
lent ISM, the effect of the hot medium pressure is essentially 
included in the turbulent velocity dispersion. 

To follow the standard technique as close as possible in 
order to facilitate comparison with previous work, the mag- 
netic and cosmic ray pressures are neglected. This is prob- 
ably the most fragile assumption. To separate variables, we 
will postpone investigating the effects of non-thermal pres- 
sure terms until section §3.5.2. 

Equation (35) must be solved under boundary condi- 
tions at high z, which are given by the properties of the 
embedding ambient medium. One can see that reasonable 
hot halo components in hydrostatic equilibrium, as that de- 
rived by Benjamin & Danly (1997), have little impact on 
the Hi scaleheight. Not only the pressure by the hot halo 
gas but also any ram pressure exerted by the intergalactic 
medium onto the disc is also thought to have an ignorable 
effect on the vertical distribution of the H I layer. 

To date, most authors have adopted these assump- 
tions to determine the distribution of the gas (e.g., Oiling 
1995, 1996; Becquaert & Combes 1997; Oiling & Merri- 
field 2000; Narayan & Jog 2002; Kalberla 2003; Narayan et 
al. 2005; Abramova & Zasov 2007). Although the limitations 
of these suppositions are well-documented in the literature, 
the reader is referred to §3.5 for a discussion. 



3.1 The method 

We have derived the vertical distribution of the H I layer fol- 
lowing the same simplifying assumptions than in the stan- 
dard procedure implemented in Oiling (1995). The system is 
assumed to be in a (perhaps turbulent) steady state which 
can be described by a Jeans-like equation. The velocity dis- 
persion tensor of the gas, which includes thermal and macro- 
scopic gas motions, is assumed to be round and constant 
(isothermal) with z as well as R. The maps of the second 
moment (or velocity dispersion) in external galaxies display 
a patchy distribution (e.g., Petric & Rupen 2007) with re- 
gions of relatively low velocity dispersion (5 km s _1 ) and 



3.2 Mass model 

In order to estimate the thickness of the H I layer, the circu- 
lar velocity should be known (see Eq. 24). Given the mea- 
sured baryon distribution, the circular velocity can be deter- 
mined. In this sense, MOND has less free parameters than 
the hypothesis of dark haloes. We consider the Galaxy to 
be composed of a bulge, a stellar disc and ISM in the form 
of Hi and H2 layers. For our Galaxy they are all based on 
well-known observations as described below. Our analysis 
treats the gas and the stars as axisymmetric distributions. 
Therefore, the resulting quantities should be considered as 
average values over azimuth. 
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3.2.1 Stellar components 

The density profile of the bulge is assumed to be described 
by a spherical Plummer-Kuzmin model: 



Pb(r) 



3M b 
4ttR?. 



1 + 



5/2 



(36) 



where the bulge mass Mb = 3.2 x 10 10 Mq and the bulge 
scalelength Rt = 2.5 kpc. Inclusion of this simple bulge re- 
produces a reasonable looking rotation curve of the Galaxy. 
The contribution of the bulge to the Hi scaleheight in the 
outer parts (R > 10 kpc) of the Galaxy is negligible. For 
instance, for a mass Mb = 1.5 x 10 10 Mq (e.g., Flynn et 
al. 2006), the scaleheight changes < 4%. 

The stellar disc is modelled as an exponential surface 
density distribution with a central surface density Eo and 
disc scalelength Rd such that: 

Y,*(R) = Y; e- R/Rd , (37) 

where S = 640.9 M pc~ 2 and R d = 3.2 kpc (Mera et 
al. 1998; Lopez-Corredoira et al. 2002; Larsen & Humphreys 
2003; Sommer-Larsen et al. 2003), implying a total stellar 
mass of 4.1 x 10 10 Mq in the disc. The local stellar sur- 
face density (at R —Rq = 8.5 kpc) turns out to be 45 M Q 
pc~ 2 , which when combined with the corresponding value 
for gas (Hi + H2) of 7 Mq pc~ 2 , gives the total local sur- 
face density inferred by most observations (e.g., Holmberg 
& Flynn 2004) . The vertical velocity dispersion of the stars 
as a function of R is taken from Lewis & Freeman (1989). 

3.2.2 Gaseous components and the flaring of the Hi layer 

The H 1 surface density (He corrected 1 ) as a function of the 
galactocentric radius is taken as an input for the disc model. 
As said before, the Hi data comes from the most recent 
comprehensive Hi survey of our Galaxy, the LAB survey 
(Kalberla et al. 2005). Here we use the average observed 
flaring (of North and South hemispheres) results for a flat 
rotation curve with v c = 220 km s _1 and Rq = 8.5 kpc (the 
current IAU standard values), to compare with the predicted 
scaleheight in MOND. 

The vertical velocity dispersion of Hi, cthi, outside the 
solar circle is hard to measure, and so not yet known. Ob- 
servations of external galaxies imply that <thi is typically 
constant along galactocentric radius with a value 7-9 km 
s _1 (e.g., Sellwood & Balbus 1999), or declines slightly with 
radius (Dib et al. 2006, and references therein). Since in the 
inner Galaxy the H I is observed to have a velocity disper- 
sion of 9.2±1 km s _1 independent of radius (Malhotra 1995) 
and since there is evidence that the velocity dispersion in the 
outer Galaxy equals the value inside the solar circle (Blitz 
& Spergel 1991), Oiling & Merrifield (2000, 2001) assumed 
a constant cthi of 9.2 km s~ as a more likely value. How- 
ever, a slightly declining velocity dispersion could also be 
admissible. In order to calculate the scaleheight we consider 
models across a range of values 7 km s _1 < cthi < 9 km s _1 , 
constant with radius. 

The H2 surface density as a function of R was taken 



1 Our analysis does include a 30% per cent contribution by mass 
from helium. 



from Wouterloot et al. (1990) and its vertical velocity dis- 
persion is 5 km s _1 (Clemens 1985). 



3.3 MOND parameters and the rotation curve 

Milgrom's formula (Eq. 1) has been used extensively to fit 
the rotation curves of external galaxies. The MOND circular 
speed v c in the disc mid-plane is related to its Newtonian 
counterpart (without any dark matter), v c ,n, by the follow- 
ing equation: 



li{x)v 2 c = v% N (38) 

(Brada & Milgrom 1995) where x = vi/(Rao) and the New- 
tonian rotation speed in the mid-plane of the Galaxy is given 
by: 
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(39) 



For a set of external galaxies with high quality rotation 
curves, Begeman et al. (1991) derived a = (1.2±0.27) x 10~ 8 
cm s~ . This has been the preferred value for MOND stud- 
ies that followed (Sanders & Verheijen 1998; Famaey & Bin- 
ney 2005). However, Bottema et al. (2002) argue based on 
Cepheid-based distance scale to UMa cluster of galaxies that 
the value of ao should be adjusted to 0.9 x 10~ 8 cm s~ 2 . 
Adopting the mass model described in §3.2 and the relation 
between the Galactic constants, O0/-R© = (27 ± 2.5) km 
s^kpc" 1 (Kerr & Lynden-Bell 1986; Reid et al. 1999), rea- 
sonable rotation curves are generated for ao in the range of 
0.8 x 10~ 8 to 1.2 x 10~ 8 cm s . Figure 1 shows the rota- 
tion curves for the interpolating functions with n = 1 and 
n — 2 and different a . The radial upper limit of 40 kpc is 
set by the reliability of H I data. At the solar point Rq = 8.5 
kpc, the circular speed lies in the range of 215 to 243 km 
s~\ conforming with the accepted range of ratio of Galactic 
constants. This implies that the model parameters assumed 
above are appropriate. Note that the differences in the rota- 
tion curves are small at large radii because the asymptotic 
value of the rotation speed is given by v c (oo) = (GMao) 1 ^ 4 
for both interpolating functions. 

We like to point out that the assumed rotation curves 
behind the theoretical and observed scaleheight curves do 
not exactly match. In the case of H I observations, the con- 
version of brightness temperature Tb , to volume density de- 
pends on the assumed distance to the Galactic centre Rq, 
the local rotation velocity vq and the outer rotation curve. 
The conventional way to go with LAB survey analysis is to 
assume a flat rotation curve of 220 km s _1 with Rq = 8.5 
kpc. The scaleheight data used in this paper is based on 
these values, whereas the rotation curves used to estimate 
the theoretical flaring are based on the adopted mass model 
of our Galaxy. Ideally, we should use rotation curves in Fig. 1 
to analyse the LAB data. However, the difference between 
the two is only 10% — 15% and we will see that, for our 
purposes here, the results are stable against such modest 
adjustments. So its effect on the scaleheight will not com- 
promise our conclusions. 



3.4 H 1 scaleheight in the Galaxy with MOND 

Using Runge-Kutta method, we solve equations (24) and 
(35) numerically as an initial value problem to obtain our 
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Figure 1. Rotation velocity curves for different ao, in units of 
10~ 8 cm s~ 2 , and for the FB (panel a) and standard (panel 
b) interpolating functions. Note that the circular speed v c is 
slightly higher for the FB interpolating function. MOND pro- 
duces a falling rotation curve in the region of interest (R > 10 
kpc). 



object of interest: density as a function of z at different R. 
Readers are referred to Narayan & Jog (2002) for a detailed 
numerical procedure. As discussed in §2.5, we can neglect £ 
as compared to 1 in Eq. (24) in the whole range of interest 
(10 kpc < R < 40 kpc) because at large radii where the L- 
term is important, £ < 0.04 (see the Appendix). In this way 
we obtain the simultaneous vertical density distributions of 
stars, Hi and H2. 

The vertical density distributions of stars and Hi for 
am = 7 km s _1 , n = 1 and ao = 1.2 x 10~ 8 cm s~ 2 arc 
shown in Fig. 2. Even though we have three disc compo- 
nents (stellar disc, Hi and H2 layers) having different ver- 
tical velocity dispersions and gravitationally coupled under 
MOND, the vertical density profiles are still very close to a 
sech 2 -function. As expected, the H2 layer is the most com- 
pact one followed by the Hi layer and the stellar disc, re- 
spectively (see Fig. 2, upper panel). The scene is different 
at large distances from the Galactic centre. The azimuthally 
averaged H2 surface density is so low there that it becomes 
incorrect to assume it as a continuous layer. It may exist as 
rare clouds and so disappear from our picture. 

By solving the equations at regular R intervals (between 
10-40 kpc), the scaleheight h, defined as the HWHM of the 
vertical Hi density distribution , can be plotted as a func- 
tion of R (Fig. 3). The Hi layer undergoes flaring due to a 
combination of decreasing total surface density and a steady 
velocity dispersion of 7 km s _1 , and therefore extends verti- 



2 Note that the definitions of zo and h are different. In this sec- 
tion, it is more convenient to use h. 
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Figure 2. Density distribution along the vertical direction for 
H I, stars and H2 in the outer parts of the Galactic disc according 
to MOND for <t H i = 7 km s -1 , a = 1.2 X 10~ 8 cm s~ 2 and 
n = 1. 



cally beyond the stellar disc. The L-term in Eq. (24) becomes 
similar to the local self-gravity term, 4tt[i~ 1 Gp, at ~ 27 kpc. 
At R = 40 kpc it turns out to be 4 times larger than the 
local self-gravity when the FB interpolating function and 
ao = 0.8 x 10~ 8 cm s~ 2 are used. While the differences in the 
flaring between the standard and the FB interpolating func- 
tions are not significant beyond R = 27 kpc, the standard 
interpolating function predicts a thicker Hi disc at 10 kpc 
< R < 20 kpc (see §3.5.3 below). For n = 2, a = 1.2 x 10~ 8 
cm s~ 2 and am = 7 km s _1 , MOND predicts a scaleheight 
of ~ 375 pc at R = 20 kpc and ~ 1625 pc at R = 40 kpc. 

Now we make a one-to-one comparison between the pre- 
dicted and the observed scaleheights of Hi from the LAB 
survey. This comparison is possible only for the H I layer be- 
cause the number counts of both stars and molecular hydro- 
gen fizzle out at such large radii. The circles in Figs. 3 and 
4 represent the H I scaleheight from the LAB survey data 
as derived in the reference model (see §3.2.1 in Kalberla et 
al. 2007). For the adopted am, the predicted MOND scale- 
height always falls below the LAB values in the entire range 
considered. Close to 10 kpc, the predicted value is almost 
60% below the observed value and at 40 kpc it is about 30% 
less. The lower ao slightly lifts the curve up but the discrep- 
ancy at the interval 10 kpc< R < 17 kpc and beyond 30 kpc 
still remains. 

The velocity dispersion in our Galaxy is, of course, un- 
certain. In the inner Galaxy, the Hi is observed to have a 
velocity dispersion of 9.2 ± 1 km s _1 , independent of radius 
(Malhotra 1995). If we assume that this value characterizes 
the turbulent motions of the gas throughout the Galaxy, the 
flaring of the H I layer can be reproduced well at R > 17 kpc 
(see Fig. 4). At 25 kpc < R < 32 kpc, the predicted flar- 
ing curve moves away from the observed one by less than 
15%. Due to the fact that the flaring in the northern hemi- 
sphere is ~ 80% larger than in the sourthern hemisphere, 
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the mentioned discrepancy of 15% can be attributed to our 
axisymmetric assumption. It is important to note that the 
L-term and the local self-gravity term vary having a different 
radial dependence and should work out in phase to explain 
the flaring curve. If the L-term were not included, MOND 
would overestimate the thickness of the H I disc beyond ~ 25 
kpc. In particular, it would be overestimated by a factor of 
5 at R — 40 kpc. The success of MOND in predicting the 
H I flaring at the outer parts is noteworthy because at such 
large radii the predictions are very robust to small adjust- 
ments on the Galactic mass model, Galactic constants, the 
value of ao or the index n of the interpolating function (see, 
for instance, Fig. 3). The thickness is essentially determined 
by the adopted cthi. If a Hi velocity dispersion of w 9 km 
s _1 is taken at face value, the observed H I data up to 40 kpc 
seems to favour the lowest acceptable values of the universal 
acceleration ao ~ 0.8 x 10~ 8 cm s~ 2 . 

In Fig. 4 we see that the observed thickness at 10 kpc 
< R < 15 kpc is still a factor of 1.7 larger as compared to the 
MOND prediction. This factor cannot be attributed to the 
approximations we made in the derivation of Eq. (24). It is 
also unlikely that systematic deviations from axisymmetry 
can account for that difference because the average thickness 
over 0° < (f> < 180° is similar to the average value over 
180° < 4> < 360° at the interval 10 kpc < R < 15 kpc 
(Levine et al. 2006a; Kalberla et al. 2007). 

It should be borne in mind that the flaring curve and 
the H I volume density as inferred from the LAB data are de- 
pendent on the Galactic constants and the adopted rotation 
curve and we pick up the v c — 220 km s _1 reference model. 
Nevertheless, we established, after some investigation, that 
the above mentioned discrepancy between the predictions 
and observations cannot be smeared out by slightly chang- 
ing the values of Rq and wq, or by adopting the rotation 
curves of Fig. 1 to analyse the LAB data. In order to sort 
through a large parameter space, we look at other param- 
eters that do not appear to be strongly constrained by the 
data; this will be done in the next section. 



3.5 Effects of incomplete modelling and other 
uncertainties 

In the previous section it was found that in the MOND 
formulation the observed disc material gives rise to too thin 
a thickness at 10 kpc < R < 15 kpc. In order to determine 
the implications of this result, we need to consider effects of 
incomplete modelling and other uncertainties. 



3.5.1 The stellar mass-to-light ratio 

As already said, the L-term is small as compared to the 
self-gravity term in the "inner" disc (R < 20 kpc) , implying 
that h oc /x/S. It is possible to have a thicker H I disc in the 
inner disc if £ is decreased simply by decreasing the stellar 
mass-to- light ratio 3 T. In fact, if T is decreased, the circular 



3 The associated error in the local column density of stars is ~ 10 
M Q pc~ 2 (Flynn & Puchs 1994; Holmbcrg & Flynn 2004). In ad- 
dition, strong departures from axisymmetry could mean that the 
local surface density may not be representative of the azimuthally 
averaged value. 



velocity and thus /i lower but the drop in jj, is not enough to 
compensate the reduction on E. In particular, if we take the 
FB interpolating function with ao = 0.8 x 10~ 8 cm s~ 2 and 
reduce T to a value such as that vq = 200 km s _1 , which 
is quite low but still compatible with the terminal velocity 
data (Famaey k, Binney 2005), h increases about 27% at 
R — 15 kpc and about 10% at R > 35 kpc. Therefore, the 
fit is globally better in the latter case but the overcorrection 
is not eliminated. 



3.5.2 Non-thermal pressure terms 

A velocity dispersion ctht of ~ 9 km s _1 may be attributed 
to the gas random motions. Potentially, there may be other 
forces helping to support the H I layer as cosmic rays and 
magnetic fields, which lead to a larger effective velocity dis- 
persion a m (e.g., Parker 1966; Spitzer 1978; Boulares & Cox 
1990). The variation with radius in the pressure support 
from magnetic fields and cosmic rays is essentially uncon- 
strained. The simplest assumption is that, like the turbu- 
lent pressure, these two terms do not vary with radius and 
so the effective velocity dispersion is a constant. Following 
Spitzer (1978), denote by as the ratio between the mag- 
netic pressure arising from the regular field Pb and the ki- 
netic pressure P g . Similarly, we introduce ab = PbjPg and 
acR = Pcn/Pg, where Pb is the pressure arising from the 
random field and Pen the pressure due to cosmic rays. The 
scaleheights of the non-thermal pressure terms are Kb, Kb 
and kcr times larger than the gaseous scaleheight. With 
this notation, the relative increase in the velocity dispersion 
is 



frm _ L | Qb | ab acR 

CHI V K B Kb KCR 



(40) 



(e.g., Oiling & Merrifield 2001). Models and observations of 
the ISM in the solar vicinity and external galaxies are com- 
patible with approximate equipartition between the energy 
density of the turbulent magnetic field and the turbulent en- 
ergy density of the gas (e.g., Zweibel & McKee 1995; Beck et 
al. 1996; Beck 2001; Fletcher & Shukurov 2001; Beck 2007), 
hence a b ~ 0.3-0.5 (see also Kalberla & Kerp 1998). Ob- 
servational constraints on the magnetic field near the Sun 
yield 1 < b/B < 3 (e.g., Fletcher & Shukurov 2001 and 
references therein), thus qb = 0.05-0.5. Equipartition ar- 
guments suggest the cosmic ray pressure to be equal to the 
total magnetic pressure, implying ocr ^ 0.5-0.7. The ver- 
tical distribution of the non-thermal components are even 
more difficult to obtain than is that of the gas components 
of the ISM. If the turbulent magnetic field has the same 
origin than the turbulent motions, Kb = 1. For the regular 
magnetic atmosphere Kb — 5-10 are reasonable values. The 
possible values of kcr can be derived by noting that the 
cosmic ray pressure is equal to the total magnetic pressure. 

Using Eq. (40) with the above range of parameters 
we obtain a m /am = 1.2—1.4. This is in accordance to 
the results of Kalberla & Kerp (1998) who conclude that 
a' m /am « 1.15 due to the turbulent magnetic field com- 
ponent, when the pressure contribution from cosmic rays is 
not included. Ratios a m /am ~ 1.4 are commonly used in 
the literature (e.g., Sellwood & Balbus 1999; Elmegreen & 
Hunter 2000). 
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Since the discrepancy in the scaleheight is 25%-70%, 
depending on the adopted values of n and T, values of 
Chi/chi = 1.12—1.3 are required to account for the thick- 
ness of the disc at 10 kpc < R < 15 kpc. We conclude 
that, given the uncertainties, simplified estimates indicate 
that non-thermal pressure support could be sufficient. Be- 
yond the optical radius ~ 15 kpc, where star formation is 
almost nonexistent, energy input into cosmic rays and mag- 
netic fields from stellar sources is likely to be unimportant. 
Hence, a decreasing cr^j with R is expected. 



3.5.3 Form of the interpolating function 

In section 2.1, it was said that analyses of the rotation 
curves of our Galaxy and other external galaxies constrain 
the exponent n of the interpolating function in the range 
0.85 < n < 2. Whilst at the far outer disc R > 27 kpc the 
thickness of the H I disc remains almost unchanged when n 
is varied between 0.85 and 2, the flaring of the disc in the 
inner disc (10 kpc < R < 20 kpc) scales approximately as 
h oc n and so that it depends on the adopted value of n. 
Once the baryonic mass distribution is known, the ratio be- 
tween scaleheights derived with n = 2 and n = 1 in the 
inner disc is: 



h(n = 2) 
h(n = 1) 



V2 



(yr+v-i) 
yi + 4^-1 



1/2 



(41) 



where y = (v c N /Rao)~ l . In the range of interest we have 
y > 1 implying that the scaleheight is about 25% larger 
when the standard interpolating function is used than in 
the FB case. Consequently, MOND requires less non-thermal 
support, a' m /am — 1-18 for our standard set of parameters 
described in §3.2, to explain the observed flaring if the stan- 
dard interpolating function (n = 2) is adopted. 



3.5.4 Outer spiral arms and warps 

It is possible that the H I layer is slightly out of hydrostatic 
equilibrium because of the warp and the presence of outer 
spiral arms. The scattering of stars and compact clouds by 
the spiral arms hardly increases their vertical velocity dis- 
persion (e.g., Jenkins & Binney 1990). However, if the Hi 
layer behaves as a smooth disc of gas, it develops a com- 
bination of a shock and a hydraulic jump; the gas may 
shoot up to higher z (e.g., Gomez & Cox 2002, 2004). Nev- 
ertheless, the arms outside the corotation radius are likely 
weak and these vertical complications are expected to oc- 
cur along a thin (a few hundred parsecs) region downstream 
the shock/jump of the arm. The observed correlation be- 
tween the maps of H I thickness and surface density -regions 
of higher surface density having a more reduced thickness- 
(Levine et al. 2006b) , suggests that hydrostatic equilibrium 
is a good approximation. By averaging over azimuth, the al- 
teration of the vertical structure by this jump/shock should 
not be significant. 

It is well-known that the H I disc of the Galaxy is 
warped. A warp is essentially a vertical m — 1 mode that is 
maintained by some external driving force. The S-shape of 
the warp does not affect appreciably the local self-gravity in 
the disc because the wavenumber of the warp k <C 2-K/h. 
Since the orbit of a parcel of the disc is inclined to its 



galaxy's equatorial plane, its height z above this plane os- 
cillates around z — with amplitude h w (R). It is important 
to know at what extent this oscillatory motion may gener- 
ate vertical compressions or decompressions and produce an 
azimuthal dependence on the thickness of the H I disc. First, 
compressions should induce a pressure gradient in the verti- 
cal direction of the disc (probably a shock wave or compres- 
sion front) and, consequently, a characteristic asymmetric 
vertical morphology is expected (Sanchez-Salcedo 2004). By 
contrast, the observed fairly symmetric distribution of the 
Galactic warp suggests that the gas layer moves vertically in 
a solid-like fashion, as if it were incompressible. This solid- 
like behaviour is expected either if the bending modes are 
long-lived eigen-modes or they are excited by gravitational 
torques 4 . Just for illustration, let us consider the tidal ex- 
citation as the mechanism responsible for maintaining the 
warp with amplitude h w (R). The gravitational perturber 
(e.g., the Magellanic Clouds) may exert a tidal stretching 
on the disc, but this effect is very small across the thick- 
ness of the Hi disc and hence we only need to study the 
restoring gravitational force that pulls the gas back to the 
disc midplane. For low amplitudes of the warp h w /R <ti 1, 
the underlying gravitational restoring force due to the disc 
is described by an harmonic potential, $(z) = u 2 z 2 /2 with 
v the vertical frequency derived in §2.5. Under this poten- 
tial, one can show that the oscillatory mode in the vertical 
direction occurs at constant thickness. In fact, denoting by 
Peq(z) the hydrostatic density profile in the absence of the 
warp, the density p(z,t) — p eq (x) with x = z — h w cos(^i + £) 
describes an exact solution of the isothermal Jeans vertical 
equation for any value of the gas velocity dispersion. Since 
this 'rigid' mode is reached in a dynamical timescale 
a net departure from Jeans hydrostatic equilibrium because 
of the warping should be small. This 'incompressible' be- 
haviour has been confirmed numerically even in the case 
that the disc evolves only through its own gravity (see, e.g., 
§4 in Masset 1997). 



3. 5. 5 Ram pressure and cold gas accretion 

The ram pressure due to infalling diffuse material from 
the intergalactic medium could compress the disc vertically 
(Sanchez-Salcedo 2006). The dynamical collapse time for the 
H 1 layer is so short that Eq. (35) is still valid but should be 
solved with appropriate boundary conditions. What we have 
found is that MOND does not require any hypothetical and 
unseen compression by smooth intergalactic accretion flows. 

The Galactic disc is probably subject to a continu- 
ous process of circulation of gas and accretion of material 
through condensing halo gas. The continuous infall of high 
velocity clouds onto the H I disc may produce enough agita- 
tion in the interstellar medium to explain the observed levels 
of turbulence at the outer parts (Santillan et al. 2007; Booth 
& Theuns 2007). Even in this case, hydrostatic equilibrium 
of the thin H I layer is still a valid premise, provided that 
the value of the turbulent velocity dispersion is taken con- 



4 Gravitational torques may be produced by a companion galaxy 
(e.g., Weinberg & Blitz 2006) or by cosmic infall accretion of 
protogalactic material in the outer halo (Ostriker & Binney 1989). 
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sistently. The scarce impact of massive high velocity clouds 
may break down hydrostatic equilibrium but only locally. 



3.5.6 Uncertainties in the surface density of the outer disc 

Our mass model does not include stellar streams as those 
detected recently (Newberg et al. 2002; Ibata et al. 2003; 
Yanny et al. 2003; Martin et al. 2005; Martmez-Delgado et 
al. 2005). Since a ring- like mass concentration produces an 
enhanced vertical acceleration in MOND, the presence of 
such structures are sufficient for explaining the relatively 
low flaring as compared to the predicted values at R ~ 27 
kpc in Fig. 4. We refrain to find the best-fitting model until 
the distribution of the mass in the stellar disc is measured 
more accurately. 



3.6 MOND versus CDM 

One is tempted to ask: which competitive model, MOND 
or dark matter, reproduces the flaring of our Galaxy with 
the fewest ad-hoc assumptions? Narayan et al. (2005) and 
Kalberla et al. (2007) have explored the shape and density 
profile of the dark matter halo using the flaring technique. 
A nearly spherical dark halo as derived from Sgr dwarf con- 
straints provides a very good fit to the observations as long 
as the dark matter density decays very fast, as r -4 (Narayan 
et al. 2005) . In this situation, the total halo mass is 2.8 x 10 
Mq, too low to be acceptable. Accretion of homogeneous, 
diffuse intergalactic gas leads to more thin discs (Sanchez- 
Salcedo 2006), thereby going in the wrong direction to ac- 
count for the observed flaring. A complex dark matter distri- 
bution (halo + disc + ring) as seen in Kalberla et al. (2007) 
may explain the thickness of the H I disc but it remains un- 
clear if this model is compatible with the orbits of Sgr and 
its stellar debris. Kalberla et al. (2007) conclude that, if 
the dark matter particles only populates a rather spherical 
component with a flattening 0.9 < q < 1.25 (Helmi 2004; 
Johnston et al. 2005; Law et al. 2005; Belokurov et al. 2006; 
Fellhauer et al. 2006), then the observed H I thickness is dif- 
ficult to explain unless the effective velocity dispersion of 
the gas increases with galactocentric distance. 



4 DISCUSSION AND CONCLUSIONS 

Observations of spiral galaxies strongly support a one-to- 
one analytical relation between the inferred gravity of dark 
matter at any radius and the enclosed baryonic mass. This 
correlation of gravity with baryonic mass can be interpreted 
as a modification of gravity. If this is not the case, then the 
success of MOND is telling us something phenomenologi- 
cal about the nature of dark matter. Besides the shape of 
rotation curves of spirals, a wide range of gravitational ef- 
fects appear in galactic dynamics, which in principle, offer 
independent restrictions on any proposed modified gravity 
scheme. Some examples include the dynamics of stars in 
the solar neigbourhood (Kuijken & Gilmore 1987), the el- 
lipticity of the X-ray emission in elliptical galaxies (Buote & 
Canizares 1994), the gravitational stability of galactic discs 
(e.g., Sanchez-Salcedo & Hidalgo-Gamez 1999), the relax- 
ation processes in stellar systems (Ciotti & Binney 2004; 
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Figure 3. Hi scalcheight vs R for the FB (upper panel) and the 
standard (lower panel) interpolating functions and am = 7 km 
s _1 . The two curves are generated for two values of ao- Lower the 
value of ao better is the fit to the observed data. This is under- 
standable since lower ao corresponds to lower values of v c which 
means less gravitational pull. Beyond 27 kpc it really deviates 
from the observation for the above model parameters, ao is in 
units of 10 -8 cm s~ 2 . 
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Figure 4. H I scaleheight vs R for the FB interpolating function 
with ao = 0.8 X 10~ 8 cm s~ 2 and a = 9 km s -1 . For comparison, 
the scalcheight for ao = 1.2 X 10~ 8 cm s~ 2 and a = 7 km s _1 is 
also shown. 
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Sanchez-Salcedo et al. 2006) or the morphology of the Galac- 
tic warp (Weinberg & Blitz 2006). In this paper a new at- 
tempt to test the self-consistency of MOND has been car- 
ried out. Since in MOND all the gravity comes from the 
disc whereas the potential in CDM models is primarily from 
a spheroidal distribution, galactic discs are expected to be 
more flattened in MOND. For the same reason, if MOND is 
true, Newtonists should only infer oblate dark haloes. In this 
sense, a clear-cut inference of a prolate dark halo would be 
enough to rule out 'classical' MOND. Interestingly, studies 
on the shape and flattenings of dark haloes can be immedi- 
ately used as tests to modified gravities. 

Our knowledge about the shape of the Galactic dark 
halo is growing rapidly. Recent studies of the tidal debris 
from the Saggitarius dwarf and its bifurcation constrain the 
flattening of the dark halo potential to be close to spher- 
ical (0.9 < q < 1.25). Contrary to the naive expectation 
that MOND should produce far more precession because all 
gravity comes from the disc, Read & Moore (2005) found 
the amount of orbital precession is nearly identical to that 
which occurs in a CDM Galactic halo with q — 0.9. 

The flaring of H I have also proven to be useful to con- 
strain the flattening of dark haloes and is extremely sensi- 
tive to the dark matter distribution and to the law of grav- 
ity; MOND is expected to squeeze the disc into a thinner 
distribution because of the enhanced vertical acceleration. 
Kuijken & Gilmore (1989) test MOND in the context of the 
Oort discrepancy and find that it tends to overcorrect some- 
what. However, they used a value of ao nearly 4 times larger 
than currently measured. Using the updated value for ao, 
McGaugh & de Blok (1998) and Famaey & Binney (2005) 
found the dynamics of stars to be consistent with MOND 
within the uncertainties. In this paper we first presented a 
framework to estimate the flaring of Hi layers in MOND. 

The MOND field equation is highly non-linear and much 
more complex that the classical Poisson field equation. Out 
of the plane it is not clear that the algebraic relation (Eq. 1) 
applies. In addition to the local self-gravity of the disc, 
MOND predicts a new term, the L-term, which has been 
shown to be equivalent to have a fictitious spherical dark 
matter halo. The self-gravity term and the L-term have very 
different radial dependences; the L-term is dominant beyond 
a certain critical radius. Both terms should work at phase to 
explain the observed flaring of galactic discs. If the L-term 
is ignored, the Galactic flaring is overestimated by a factor 
5 at the last measured point R = 40 kpc. 

Using this framework, an axisymmetric model was used 
to derive the mean Hi flaring of the Milky Way predicted 
in MOND. Our Galaxy is particularly interesting because 
we have information about some parameters, such as the 
mass-to-light ratio, that are very difficult to obtain in other 
galaxies. It is by no means obvious that given the measured 
baryonic mass of the Milky Way with essentially no free 
parameters, MOND is able to explain the flaring of the H I 
disc in an extended radial range, between 10 and 40 kpc. 
Note, for instance, that so far the Galactic flaring is only 
understood in a dark matter scenario if the Galaxy contains 
a massive dark disc with a ring (Kalberla et al. 2007; de 
Boer et al. 2007). 

In order to compare the observed gas-layer flaring to the 
predictions, the surface density of the baryonic components 
(stars and gas) and the gas velocity dispersion need to be 



known. Since the baryonic distribution is well constrained 
from observations, the procedure is essentially with one free 
parameter am between 7 km s _1 and 9 km s _1 independent 
of radius. The results depend on the choice of the Galactic 
constants, since these values affect the gas distribution as 
inferred from observations. Following the IAU recommenda- 
tions we take Rq — 8.5 kpc and Vq = 220 km s" 1 and we 
have verified that, given the uncertainties, the results are 
robust to slight changes on the values of Rq and Vq . 

For our reference set of parameters (local surface den- 
sity of stars of 45 M Q pc -2 and Rd = 3.2 kpc) , we have found 
that MOND successfully reproduces the most recent and ac- 
curate flaring curve beyond R = 17 kpc for both standard 
and FB interpolating functions, provided the vertical veloc- 
ity dispersion of Hi is taken ~ 9 km s _1 as that observed in 
the inner Galaxy (Malhotra 1995; see also Blitz & Spergel 
1991). In the interval 10 kpc < R < 15 kpc, the observed 
scaleheight is 40% and 70% in excess of that predicted by 
MOND for the standard (n = 2) and the FB (n = 1) inter- 
polating functions, respectively. Clearly, the reduced thick- 
ness predicted by MOND is unsatisfactory. We have verified 
whether uncertainties in the assumed parameters can am- 
milorate this discrepancy without altering the good fit at 
large radii (R > 17 kpc). Interestingly, MOND fails at those 
radii where the dominant term in the field equation is not 
the confining force of the L-term but the local self-gravity 
component, thus h oc /i<r 2 /E. Therefore, a reduction on the 
adopted local stellar surface density or, equivalently, on T 
would result in a thicker H I disc in the region of interest. We 
have seen that by adopting the smallest T compatible with 
observations, h increases about 27% in the inner disc. On the 
other hand, the scaleheight also depends on the orbital ac- 
celeration vl/R through /n. A reduction of the adopted value 
for Rq goes in the good direction to mitigate the overcor- 
rection but the effect is rather irrelevant. For instance, if we 
assume Rq — 8 kpc instead of 8.5 kpc, h varies only ~ 3% 
at 10 kpc < R < 15 kpc. 

The simplest way to reconcile predictions with observa- 
tions is to argue that ordered and small-scale magnetic fields 
and cosmic rays could contribute to support the disc within 
the optical disc (R < 15 kpc). Depending on n and T, the 
effective velocity dispersion required to match observations 
should be a factor ~ 1.12-1.3 larger than the turbulent ve- 
locity dispersion. We have seen that, in principle, magnetic 
fields and cosmic rays could account for this. An enhanced 
effective velocity dispersion within the optical disc (R < 15 
kpc) is reasonable in a scenario where supernovae explosions 
feed up turbulence. An enhancement a m /am = 1-3 required 
to explain the flaring when n = 1 and a local stellar density 
of = 45 Mq pc~ is likely on the high side. Therefore, 
we can say that MOND satisfactorily explains the flaring 
provided that n > 1 and E* < 45 M pc~ 2 . 

Very little further progress is likely to be made without 
further empirical constraints on the effects of non-thermal 
terms. If we learn in future that the non-thermal pressure 
forces from cosmic rays and magnetic fields have no role to 
play in the vertical support of H I discs (e.g., if the filling fac- 
tor of H I is very low), then studies of the vertical distribution 
of cold gas in our Galaxy could either rule out MOND or, 
if MOND holds true, pose strong constraints on the stellar 
mass-to-light ratio of the Galactic disc and the form of the 
interpolating function. So far, this analysis is premature be- 
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cause we believe that a certain level of non-thermal support 
is more realistic. 

The maps of H I surface density and scaleheight in our 
Galaxy are far from being axisymmetric. Nevertheless, the 
narrow anticorrelation between Hi surface density and the 
thickness of the H I layer (see Figs. 1 and 4 in Levine et 
al. 2006b) is expected when the self-gravity of the gas is sig- 
nificant in negotiating the hydrostatic equilibrium, as occurs 
naturally in MOND. Given the present uncertainties in the 
vertical support by magnetic fields and cosmic rays, we con- 
clude that MOND can plausibly explain the Galactic flaring 
and cannot be excluded as an alternative to dark matter. 
In order to discriminate between MOND and dark matter, 
we plan to extend the analysis for a sample of galaxies of 
different types, including dark dominated LSB and dwarf 
galaxies, with accurate determinations of the rotation curve 
and H I flaring. 
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APPENDIX A: ESTIMATING THE 
MAGNITUDES OF AND £ M 

In an isothermal self-gravitating disc with one dimensional 
dispersion a, the vertical acceleration at a scaleheight h is 
1.5a 2 /h and the radial acceleration is v 2 /R. Thus the ra- 
tio between the vertical and radial accelerations is = 
1.5(<7 /v c ) 2 (R/ h) . In our Galaxy, can be computed for 
the Hi layer to be 0.14 and 0.05 at 15 and 40 kpc, respec- 
tively. At regions where <C 1, i.e. h > 1.5(a/v c ) 2 R, the 
radial component determines the kinematics of the disc. For 
a galaxy like the Milky Way, this occurs at R > 6 kpc for 
the stellar disc and at R > 1 kpc for the Hi disc. At these 
radii, we have 
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where S(^) = f* z p(z')dz' is the surface density within 
a slab of thickness 2z, and v 2 (R,z) = R(d&/dR) is the 
squared planar circular velocity. 

In order to evaluate £ M , we need to estimate dg/dz at 
z ~ h. It is simple to show that 
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where we have used Eq. (Al) and the relation a 2 
■KpT 1 GTi oa h. On the other hand 
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Combining the above equations, we obtain: 
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For illustration, we may compute £ for the Milky Way at 
R = 15 kpc and R = 40 kpc. At R = 15 kpc, £ ^ 0.15 and 
decreases to £ £s 0.0016 at 40 kpc. 



